#ifndef __MATRIZ__
#define __MATRIZ__
#include <cassert>
#include <iostream>
#include <cstdlib>
#include <math.h>
#include <fstream>
#include "./matrix.h"

using namespace std;
using namespace math;

//#define delta 0.1
class Matriz {
    public:

    /**
     * Constructores
     */
    Matriz();
    Matriz(const unsigned int);
    Matriz(const Matriz&);
    ~Matriz();

    /**
     *operadores
     */

    /** deja en este parametro la suma de este con el otro */
    void suma(const Matriz&);
    /** deja en este parametro la resta de este con el otro */
    void resta(const Matriz&);
    /** devuelve el vector de largo dim del resultado de la multiplicacion*/
    void producto(const Matriz&);
    /** deja en esta matriz el resultado del producto por un escalar*/
    void producto(const double);
    /** Pongo los nombres de las var xq es mas delcarativo en este caso
        lo que hace es calcualr Ax=b modificando x y devuelve true si pudo */
    bool resolver(const double* b, double* x) const;
    /** Recibe un vector donde va a dejar los autovalores correspondiente a 
    	cada autovector devuelto en la Matriz resultado */
    Matriz autocalcular(double*,double) const;
    /** Calcula householder de esta matriz, dando como resultado la Q, y dejando
        en el parametro la R */
    Matriz householder();
    /** transpone esta matriz */
    void transponer();
    /** norma 2 de un vector */
    double static norma2(double * vector, int largo);

    /* La inversa de la matriz */
    void inversa();

    /** calcula la transformacion caracteriztica de cierta matriz **/
    void transformacionCaracteristica(const char* entrada, double* salida, int cantidad) const;

    friend Matriz operator+(const Matriz&, const Matriz&);
    friend Matriz operator-(const Matriz&, const Matriz&);
    friend Matriz operator*(const Matriz&, const Matriz&);
    friend Matriz operator*(const double&, const Matriz&);
    friend Matriz operator*(const Matriz&, const double&);

    /**
     * getters
     */

    /** devuelve un id unico corresdiente a la posicion del elemento en la matriz */
    double operator()(unsigned int, unsigned int) const;
    double operator()(unsigned int) const;
    unsigned int size() const;
    /** devuelve en el puntero del parametro el vector con la diagonal */
    void getDiagonal(double * diag) const;

    /**
     * setters
     */
    void setPos(unsigned int, unsigned int, double);
    /** transforma esta matriz en una identidad */
    void setIdentidad();

    /**
     * comparaciones y asignacion
     */
    bool operator==(const Matriz&) const;
    bool operator!=(const Matriz&) const;
    Matriz& operator=(const Matriz&);
    void load(const char*); 
    void store(const char*);
    private:
    double sumaInferior() const;
    
    /** obtiene el paso i de la matriz Q en la factorizacion QR */
    Matriz getQi(int) const;
    unsigned int dim;
    double * matriz;

    void swapRows(unsigned int i, unsigned int j);
    void mostrar() const;
};

ostream& operator<<(ostream& out, const Matriz& m) {
    out << "[";
    for(unsigned int i=0;i<m.size();i++) {
        for(unsigned int j=0;j<m.size();j++) {
            out << m(i,j) ;
            (j<m.size()-1)?out<< '\t':out;
        }
        
        (i<m.size()-1)?out<<";":out;
    }
    out << "]";
    return out;
}

void Matriz::mostrar() const {
    cout << endl;
    for(unsigned int i=0;i<this->size();i++) {
        for(unsigned int j=0;j<this->size();j++) {
            cout << (*this)(i,j) ;
            cout<< '\t';
        }
        cout<<endl;
    }
}

void Matriz::load(const char* filename){
	ifstream input(filename, ios::in);
	double buffer;
	assert(input.is_open()); //chequeo que el archivo haya sido abierto correctamente
	for(int i = 0; i < this->dim; i++){
		for(int j = 0; j < this->dim; j++){
			input >> buffer;
			this->setPos(i,j,buffer);
		}
	}

}


void Matriz::store(const char* filename) {
	
	fstream file(filename, ios::out);	
	assert(file.is_open()); //chequeo que el archivo de salida haya sido bien abierto
	for(unsigned int i=0;i<this->size();i++) {
		for(unsigned int j=0;j<this->size();j++) {
			file << (*this)(i,j) ;
			file << ' ';
		}
		file<<endl;
	}
	file.close();
}
/***********************************
 * CONSTRUCTORES ******************
 *********************************/

Matriz::Matriz(){
    this->matriz = NULL;
    this->dim = 0;
}

Matriz::Matriz(unsigned int dim){
    assert(dim > 0);
    this->dim = dim;
    unsigned int tamanio = dim*dim;
    this->matriz = new double[tamanio];
    //inicializar en 0.
    for(unsigned int i=0;i<tamanio;i++)
    	this->matriz[i] = 0;

}

Matriz::Matriz(const Matriz& otra){
    if(this != &otra){
        this->dim = otra.dim;
        this->matriz = new double[this->dim*this->dim];
        for(unsigned int x = 0; x < this->dim*this->dim; x++){
            this->matriz[x] = otra.matriz[x];
        }
    }
}

Matriz::~Matriz(){
    if(this->matriz != NULL)
        delete [] this->matriz;
    this->matriz = NULL;
}

/***********************************
* OPERADORES         ***************
************************************/

/** deja en este parametro la suma de este con el otro */
void Matriz::suma(const Matriz& otra){
    assert(size()==otra.size());
    for(unsigned int i=0;i<size();i++)
    	for(unsigned int j=0;j<size();j++)
    		this->setPos(i,j, (*this)(i,j)+otra(i,j) );
}

/** deja en este parametro la resta de este con el otro */
void Matriz::resta(const Matriz& otra){
    assert(this->size()==otra.size());
    for(unsigned int i=0;i<this->size();i++)
    	for(unsigned int j=0;j<this->size();j++)
    		this->setPos(i,j, (*this)(i,j)-otra(i,j));
}

/** devuelve el vector de largo dim del resultado de la multiplicacion*/
void Matriz::producto(const Matriz& otra){
    assert(this->size()==otra.size());
    for(unsigned int i=0;i<size();i++) {
        //guardo la fila nueva tmpmente para no sobreescribri nada
    	double aux[size()];
    	for(unsigned int j=0;j<size();j++) {
	    	double val = 0;
	    	for(unsigned int k=0;k<size();k++)
	    		val += ( (*this)(i,k) * otra(k,j) );
    		aux[j] = val;
    	}
    	for(unsigned int j=0;j<size();j++)
    		this->setPos(i,j,aux[j]);
    }
}

/** deja en esta matriz el resultado del producto por un escalar*/
void Matriz::producto(const double lambda){
    for(unsigned int i=0;i<this->size();i++)
    	for(unsigned int j=0;j<this->size();j++)
    		this->setPos(i,j, (*this)(i,j)*lambda);
}

bool Matriz::resolver(const double* be, double* x) const{
    //se copian los valores para modificarse
    Matriz aux = (*this);
    double b[this->size()];
    for(unsigned int i=0;i<this->size();i++)
        b[i]=be[i];
    unsigned int i=0,j=0;
    //1)Triangulamos la matriz junto con b
    while(i<aux.size() && j<aux.size()) {
        //1.1)Hacemos el pivoteo parcial
        unsigned int maxRow = i;
        for(unsigned int k=i+1;k<aux.size();k++) {
            if(fabs(aux(k,j))>fabs(aux(maxRow,j))) {
                maxRow=k;
            }
        }
        if(aux(maxRow,j)==0) return false;
        //1.2)Cambiamos la fila por la de maximo valor en la columna actual
        aux.swapRows(i,maxRow);
        double tmp=b[i];
        b[i]=b[maxRow];
        b[maxRow]=tmp;
        //1.3)Hacemos el siguiente paso de la triangulacion
        for(unsigned int k=i+1;k<aux.size();k++) {
            double valAux=aux(k,j)/aux(i,j);
            for(unsigned int l=j;l<aux.size();l++) {
                if(l==j) {
                    aux.setPos(k,l,0);
                } else
                    aux.setPos(k,l,aux(k,l)-valAux*aux(i,l));
            }
            b[k]=b[k]-valAux*b[i]; //tmb en b!
        }
        i++;
        j++;
    }
    //Resolvemos la matriz ya triangulada
    for(int k=aux.size()-1;k>=0;--k) {
        double val=b[k];
        for(unsigned int l=k+1;l<aux.size();l++) {
            val-= ( aux(k,l)*x[l] );
        }
        x[k]=val/aux(k,k);
    }
    return true;
}

/* Obtiene los autovalores y autovectores de la matriz
 * utilizando el metodo QR. Deja en autovalores, los autovalores de la matriz
 * y devuelve la Matriz que representa la basde de autovectores corresponeinde
 * en el order de autovalores
 */
Matriz Matriz::autocalcular(double* autovalores, double delta=0.1) const{
    Matriz A = (*this);
    double deltaTemp = delta;
    //Iteramos mientras A no sea casi triangular superior
    Matriz autovectores(this->size());
    autovectores.setIdentidad();
    while(A.sumaInferior() > deltaTemp) {
        //SIENDO A=Q*R dejamos en Q la matriz ortogonal y en A la matriz triangular superior(R)
        Matriz Q = A.householder();
        autovectores = autovectores*Q;
        //Y APLICAMOS EL METODO A=R*Q
        A = A*Q;
    }
    //Dejamos en autovalores la diagonal de A, que se corresponde
    //con sus autovalores, en orden decreciente segun el absoluto de los mismos
    A.getDiagonal(autovalores);
    return autovectores;
}

/* Devuelve la suma de los elementos de los absolutos
 * de todos los elementos que esten por debajo de la diagonal de A */
double Matriz::sumaInferior() const {
    double val=0;
    for(unsigned int i=1;i<size();i++) {
        for(unsigned int j=0;j<i;j++) {
            val+= fabs((*this)(i,j));
        }
    }
    return val;
}

/* Aplica el metodo de householder a la matriz actual "A"
 * Deja en A la matriz triangular superior R y devuelve la
 * matriz ortogonal Q talque A=Q*R 
 */
Matriz Matriz::householder(){
    /*Para obtener el metodo lo que se hace es
      ir reflejando la matriz para eliminar cada columna
      utilizando las matrices de householder. Hacemos esto
      hasta obtener Qn*..*Q1*A = R. Esa es la R de la factorizacion,
      y a Q la obtenemos haciendo Q^t = Qn*..*Q1 */
    Matriz Q(this->size());
    Q.setIdentidad();
    for(unsigned int i=0;i<size()-1;i++) {
        Matriz Qi = this->getQi(i); //La matriz que pone en 0 la i-esima columna
        Q = Qi*Q;//Q_ = Qi*..*Q1
        (*this) = Qi*(*this); //THIS=R=Qi*..*Q1*A=Qi*R_{i-1}
        //elimino errores numericos
        //Dejo en 0 toda la i-esiama columna (salvo el primer elemento)
        for(unsigned int j=i+1;j<size();j++)
            this->setPos(j,i,0);
    }
    //Y por ultimo transponemos Q=Qn*..*Q1
    Q.transponer();
    return Q;
}

/* Pasa de la matriz A a la matriz A^t */
void Matriz::transponer() {
    for(unsigned int i=0;i<size();i++) {
        for(unsigned int j=0;j<i;j++) {
            double temp = (*this)(i,j);
            this->setPos(i,j, (*this)(j,i));
            this->setPos(j,i,temp);
        }
    }
}

/* Devuelve la i-esima Matriz de Householder, para poner en 0
 * los elementos de la i-esima columna debajo de la i-esima fila */
Matriz Matriz::getQi(int columna) const {
    /* El metodo consite en reflexiones para obtener una matriz
       ortogonal talque esa matriz por la Matriz A, deje a A con 0
       en la columna deseada. Para realizar esto lo que hacemos es,
       dado la columna y fila i, nos quedamos con el vector que 
       que va desde [aii,..,ain]. a este vector le restamos en su
       primer componente, la norma del vector. A este nuevo vector
       lo llamamos u=[aii/norm,ai(i+1),..,ain]^t..
       A partir de aqui para obtener la matriz de householder, lo que hacemos
       es realizar la operacion Qi = I-2/norm(u)^2*u*u^t, y extenedmos la matriz
       para dejarla de tamno de la orignial, con 1 en la diagonal de los nuevos lugares*/
    int tamano = size()-columna; //el tamano del vector u
    double vectorResultante[tamano]; //este va a ser u
    double normaVector=0;//esta va a ser la norma del subvector [aii,..,ain]
    bool seguir = true;
    for(int i=0;i<tamano;i++) {
        vectorResultante[i] = (*this)(columna+i,columna);
        normaVector += vectorResultante[i]*vectorResultante[i];
    }
    //esto es FUNDAMENTAL, aca es donde se hacia las cosas con 0 que quedaban horrible
    //cuando la norma vector = al cuadrado del primer elemento del vector(ie, no hay q reflejar nada)
    //entoces en este caso se devuelvce la matriz Identidad! Veo dsps como resolver esto del 0.00001
    if((normaVector)-(vectorResultante[0]*vectorResultante[0]) < 0.0001)  {
        Matriz ret(size());
        ret.setIdentidad();
        return ret;
    }
    //aca vectorResultante = [aii,..,ain]
    normaVector = sqrt(normaVector); //normaVecotr = ||[aii,..,ain]_2]||
    vectorResultante[0] = vectorResultante[0] - normaVector; //vectorResultante = u
    //ahora calculamos la norma2 al cuadrado
    //de este vector y la dejamos en normaVector
    normaVector = 0;
    for(int i=0;i<tamano;i++) {
        normaVector += vectorResultante[i]*vectorResultante[i];
    }
    //Creamos Qi que es la matriz a devolver, esta tiene
    //entre de 0..columna-1:1 en la diagonal, 0 furera
    //entre columna..size() la submatriz I-2/norm(u)^2*u*u^t
    Matriz Qi(size());
    Qi.setIdentidad();//tenmos la I, ahora falta restar
    for(int i=0;i<tamano;i++) {
        for(int j=0;j<tamano;j++) {
            double restar;
            restar = vectorResultante[i]*vectorResultante[j];//esto es u*u^t en la pos adecuada
            restar = restar*2/normaVector;
            Qi.setPos(columna+i,columna+j, Qi(columna+i,columna+j)-restar); //aca queda para cada punto
        }
    }
    //y listo
    return Qi;
}



Matriz operator+(const Matriz& a, const Matriz& b){
    assert(a.size() == b.size());
    Matriz res(a);
    res.suma(b);
    return res;
}

Matriz operator-(const Matriz& a, const Matriz& b){
    assert(a.size() == b.size());
    Matriz res(a);
    res.resta(b);
    return res;
}

Matriz operator*(const double& lambda, const Matriz& a){
    Matriz res(a);
    res.producto(lambda);
    return res;
}


Matriz operator*(const Matriz& a, const double& lambda){
    return lambda*a;
}

Matriz operator*(const Matriz& a, const Matriz& b){
    assert(a.size() == b.size());
    Matriz res(a);
    res.producto(b);
    return res;
}

/* Calcula la transformacion caracteorisitica de la matriz this, 
 * con el vector entrada y devuelve en el vector salida segun cantidad
 * de autovalores a usar
 */
void Matriz::transformacionCaracteristica(const char* entrada, double* salida, int cantidad) const {
    /* La idea es para cada COLUMNA i en [0..cantidad) (autovector) de this
    la multiplicamos(t) por la entrada y dejamos el valor en la i-esima
    posicion de salida */
    double aux[size()];
    for(int i=0;i<size();i++) {
    	unsigned char valAux = (unsigned char)entrada[i];
    	aux[i] = (double)valAux;
    }
    
    for(int i=0;i<cantidad;i++) {
        double valor=0;
        for(int j=0;j<this->size();j++) {
            valor+= ((*this)(j,i)*aux[j]);
        }
        salida[i] = valor;
    }
}


double Matriz::norma2(double * vector, int largo){
    double res = 0;
    for(int i = 0; i < largo; i++){
        res += vector[i]*vector[i];
    }
    return sqrt(res);
}

void Matriz::inversa()  {
    matrix<double> image(size(),size());
    for(int i=0;i<size();i++) {
        for(int j=0;j<size();j++) {
            image(i,j) = (*this)(i,j);
        }
    }
    image = !image;
    for(int i=0;i<size();i++) {
        for(int j=0;j<size();j++) {
            this->setPos(i,j, image(i,j));
        }
    }
}

/**********************************
* GETTERS    **********************
**********************************/
/** devuelve el valor en la posicion (fila,columna) */
double Matriz::operator()(unsigned int fila,unsigned int columna) const{
    assert(fila < this->dim && columna < dim);
    return this->matriz[fila*this->dim + columna];
}

double Matriz::operator()(unsigned int id) const{
    return this->matriz[id];
}


/** devuelve el tamanio de la matriz */
unsigned int Matriz::size() const{
    return this->dim;
}

/* Pone en diag la diagonal de la matriz */
void Matriz::getDiagonal(double * diag) const {
    for(int i=0;i<this->size();i++) {
        diag[i] = (*this)(i,i);
    }
}
/**********************************
* SETTERS    **********************
**********************************/
/** setea el valor en la posicion (fila,columna) */
void Matriz::setPos(unsigned int fila, unsigned int columna, double valor) {
    assert(fila < this->dim && columna < dim);
    this->matriz[fila*this->dim + columna] = valor;
}


/* Dada una matriz, le pone 1 en toda la diagonal
 * y 0 en todo el resto
 */
void Matriz::setIdentidad(){
    for(int i=0;i<size();i++) {
        for(int j=0;j<size();j++) {
            this->setPos(i,j,0);
        }
        this->setPos(i,i,1);
    }
}

/*********************************
* COMPARACION Y ASIGNACIONES *****
*********************************/
bool Matriz::operator==(const Matriz& other) const {
    if( this != &other) {
        if(this->size() != other.size()) {
            return false;
        } else {
            bool igual = true;
            for(unsigned int i = 0; i < this->size() && igual; i++) {
                for(unsigned int j = 0; j < this->size() && igual; j++) {
                    igual = igual && (*this)(i,j) == other(i,j);
                }
            }
            return igual;
        }
    }
    return true;
}

bool Matriz::operator!=(const Matriz& other) const{
    return !((*this) == other);
}

Matriz& Matriz::operator=(const Matriz& other){
    if( this != &other) {
        this->dim = other.size();
        if( this->matriz != NULL)
            delete [] this->matriz;
        this->matriz = new double[this->dim*this->dim];
        for(unsigned int x = 0; x < this->dim*this->dim; x++)
            this->matriz[x] = other.matriz[x];
    }
    return *this;
}

void Matriz::swapRows(unsigned int i, unsigned int j) {
    assert(i<this->size() && j<this->size());
    if(i!=j) {
        for(unsigned int k=0;k<this->size();k++) {
            double tmp = (*this)(i,k);
            setPos(i,k,(*this)(j,k));
            setPos(j,k,tmp);
        }
    }
}

#endif
